New approach to Monte Carlo calculation of buckling of supercoiled DNA loops 
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The short supercoiled circular DNA molecules are shown to be glassy systems and canonical 
Metropolis Monte Carlo simulations of the systems tend to get stuck in local meteistable energy 
basins. A novel Monte Carlo algorithm is developed to alleviate the problem of "ergodicity break- 
ing" of the glassy systems, in which the Markov process is driven by an explicitly analytic weight 
factor with enhanced probability in both low- and high-energy regions. To characterize the degree 
of puckering of the supercoiled DNA loops, a new quantity of aplanarity is introduced as the short- 
est principal axis of configurational ellipsoid of DNA. With the suggested Monte Carlo method, 
the quantitative correlation between supercoiling degree and buckling of DNA is attained. With 
supercoiling stress increasing, the conformational transition from a circle to mono-, diplo- or triple 
interwound superhelical structure will take place in a successive but decreasingly abrupt mode. 
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The onset of supercoiling and supercoiling transitions 
of closed circular deoxyribonucleic acid (DNA) are in- 
triguing phenomena from both a theoretical and a biolog- 
ical viewpoint. It is generally believed that proper super- 
coiling strain is crucial for the DNA molecule to pack and 
fold into chromosomes or to interact with specific protein 
PI . In the simplest mechanical description, in which one 
neglects all molecular inhomogeneities, a double-helical 
DNA molecule can be viewed as an elastic rod, the elas- 
tic energy of which can be decomposed into bending and 
torsional parts [0. These two energy terms are usually 
decoupled since the bending energy depends on local cur- 
vature of the center axis, while torsional energy is the 
function of displacements of twist of the DNA basepairs 
from the equilibrium positions. For a torsionally relaxed 
small loop, (e.g., nicked DNA polymer with one broken 
strand), the torsional energy of DNA vanishes, and DNA 
molecule has the minimum bending energy of a flat cir- 
cle. The effective coupling between bending and torsional 
energies rises for covalently closed DNA polymers, since 
the number of times two strands of the DNA duplex are 
interwound, i.e., the linking number Lk, is of topological 
invariance. The planar circle would no longer stand as a 
stable minimum of the total elastic energy when the link 
number deficit/excess is beyond some threshold which 
is a function of bending and torsional stiffness of DNA 
polymer 0]. In fact, the spatial configuration of circu- 
lar DNA depends upon the competition between bending 
and twist energies. Planar DNA microcircles will pucker 
if the decrease in torsional energy through changing the 
twist of the basepairs exceeds the increase of bending 
energy caused by writhing of the polymer axis. 

The buckling of short supercoiled DNA loops of 1000 
basepairs (bp) was computationally studied by Schlick 
et al. with deterministic techniques of energy minimiza- 
tion and molecular dynamics (MD) 0]. A catastrophic 
buckling transition from the circle to the conformation of 
figure-8 was presented with increasing supercoiling strain 
of DNA, coinciding with earlier analytical calculations 
[pi . Metropolis Monte Carlo (MC) approach to this buck- 
ling transition of DNA circle of 468 bp was performed by 



Gebe and Schurr, through tracing the writhing of the 
molecule 0. 

In this Rapid Communication, we study the buckling 
dynamics of DNA loops of 168 bp through a new defined 
puckering degree (see below), by Monte Carlo simulation 
of discrete wormlike chain model [|6|. The deformation 
energy of DNA rod is then approximated by the har- 
monic bending and twisting components. In each update, 
an interval subchain containing arbitrary amount of links 
is rotated around the straight line connecting the vertices 
bounding the subchain. The rotating angle is randomly 
taken in a interval chosen so that about half of the pro- 
posed moves are accepted. The excluded-volume effects 
and electrostatic repulsive interaction are incorporated 
by a hard-cylinder potential with an effective diameter 
(d > 2 nm) in the way the free energy of DNA will be infi- 
nite when the distance of any non- adjacent parts of DNA 
is less than d. However, as shown in Fig. nk, the conven- 
tional Metropolis simulation |7| of short DNA polymer of 
168 bp with Boltzmann weight factor tends to get stuck 
in some configurations with local metastable minimum 
energies. 

In fact, the energy landscape of supercoiled DNA chain 
is charactered by numerous local minima separated by 
energy barriers. At length scales comparable to the 
double-helix repeat of 3.4 nm (« 10.5 basepairs) or the 
diameter of 2.0 nm H, the pairing and stacking enthalpy 
of the bases make the polymer remarkably rigid and the 
energy barriers significant high compared with Boltz- 
mann weight factor. Since the probability of a canon- 
ical Metropolis procedure to cross the energy barrier of 
height AE is proportional to exp{— AE / ksT) , the simu- 
lations therefore tend to get trapped in some local energy 
basins (Fig. |l|a), although Metropolis Monte Carlo sam- 
pling is usually more effective than molecule dynamics 
simulation for conformational changes which jumps to 
different area of phase space [|| . For the glassy systems 
such as short DNA loops, different MC and MD simula- 
tions within finite CPU time and sweeps may get stuck 
in different energy basins, which renders the calculations 
of physical quantities unreliable. One of the main aims 



of present work is at alleviating this problem of "ergodic- 
ity breaking (EB)" in Monte Carlo simulations of glassy 
systems. 



run, which efRciently keeps the simulation from getting 
stuck in local minima. 
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FIG. 1. The energy of DNA loop of 168 basepairs with 
linking number deficit ALfc = 5, as a function of computer 
time: (a), for canonical Metropolis method, the simulation is 
getting stuck in a local energy basin at first 700,000 steps and 
in another local energy minimum afterwards; (b), for present 
sampling with new weight factor of Eq. (|l|), much larger fluc- 
tuations are implemented, which keeps the simulation from 
getting trapped in local energy basins. 

Technically, the above-mentioned EB problem is due to 
the fact that Boltzmann weight factor drops too quickly 
(exponentially) with system energy E. Therefore, the 
probability of the Markov process to jump over the high 
energy barriers is exponentially damped. In the follow- 
ing, we consider a new weight factor 



w{E) = e 



^ -E+{V2/a)\E-{E)\ 



(1) 



where (E) is the averaged energy of system, and a^ 
(= np /2) is the mean squared deviation of energy of the 
canonical thermodynamic system, and np the number 
of degree of freedom. For the discrete closed wormlike- 
chain with N links, the number of degree of freedom is 
Up = 2N — 6. Here and after the energy and rigid- 
ity parameters are scaled by fc^T, so the Boltzman in- 
verse temperature parameter /3 = l/fcsT is omitted in 
our equations. 

With the introduction of factor of cxp[{y/2/(T)\E — 
{E)\], the sharp peak of energy distribution of canoni- 
cal ensemble is damped and the important sampling in 
low-energy region is reinforced. Especially, within the 
new weight factor of Eq. (|l|) , the probability in higher- 
energy region is enhanced exponentially, which makes the 
simulation escape from local energy basins easily. In Fig. 
nb we present the MC time series of energy of the same 
DNA system but with the new weight factor of Eq. (|^). 
Indeed, the simulation of the new artificial ensemble cov- 
ers a much wider energy range than that of the canonical 
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FIG. 2. The probability distributions of energies of DNA 
loops, which is calculated by reweighting the artificial ensem- 
ble with Eq. (Q) . After translating the energy to the averaged 
value (E) , the distributions of different linking number deficit 
ALk fall on the same curve. The results of the simulation 
by entropic sampling method pT|-|l4| coincide with that by 
present approach, but its CPU time is twice more than the 
latter. 

In the artificial ensemble of (|^), each configuration 
of energy E of DNA loop, in fact, represents n{E) 
(= exp[—{^/2/a)\E — {E)\]) ones of the real thermody- 
namic system. Through reweighting the artificial sample 
Pi, i.e. 



(A) = 
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where A^swoop is the number of MC sweeps of the ar- 
tificial sample, we can obtain the expectation value or 
probability distribution of the considered quantity A in 
the real physical system. As an illustrated example, we 
show in Fig. 2 the probability distributions of DNA en- 
ergies, which are calculated from the artificial sample of 
2,000,000 MC runs. As expected, after the translation of 
the DNA energies to their averaged values {E), all the 
distributions with different supercoiling strains precisely 
fall onto the same curve. In fact, this scaling behavior of 
the energy distribution of canonical ensemble is the foun- 
dation that we could use the same compensative factor 
in Eq. (mj to alleviate the EB problem in different energy 
systems. 

It should be mentioned that the EB difficulty has been 
longly standing in the Monte Carlo simulations of glassy 



systems. The most often-used technique of earher deal- 
ing with this problem in first-order phase transitions [[ll| , 
pro tein folding problem |12| and other glassy systems 
||l3| , is entropic sampling method [Ijjl^, in which the 
spectral density of the considered system is calculated 
numerically by an iterative procedure so that the simu- 
lation can be performed as a random walk in a desired 
range of phase space. Depending on the rate of conver- 
gence of iteration and the size of simulated energy range, 
the determination of weight factor is usually tedious and 
very much time-consuming p^-p^. To get an approxi- 
mate flat distribution of energy spectral in the simulation 
of our short DNA loops, for example, we need to spend 
more than 50% of the whole CUP time to achieve this 
task. In Fig. 2 we also show the probability distribution 
of DNA energies of the canonical ensemble obtained from 
the entropic sampling technique for ALfc — and 5 re- 
spectively. The results coincide with that by the present 
method. 

The only parameter to be determined in Eq. (n^) is 



the averaged energy of the system. In principle, such 
estimates can be found in an simple iterative way [[l6[ . 
One first sets an initial guess of the averaged energy 
and performs a simulation with small number of Monte 
Carlo sweeps, and gets a new value of the averaged en- 
ergy which is a better estimator for the real energy. 
One can run the simulation using the new estimator of 
(-E) to get a newer one and iterate this process until 
(E) w ^J {E'^) — 77.^/2. According to our calculations, 
the averaged energy converges to the stable value very 
quickly, and just a few times of iterations are enough to 
get a precise estimator of the averaged energy. In Ta- 
ble I is listed the averaged energies of the DNA loops of 
some different supercoiling restrains, which are obtained 
by 3 iterations. The estimation of averaged energy can 
be further facilitated by information of its ground state 
value Eq. For example, the ground state of DNA loop for 
AL/c = corresponds to the flat circle, energy of which 
is 18.14/cBr in our system. So we can directly calculate 
the averaged energy by {E) = Eq + np/2 = i^AAksT. 



TABLE I. The averaged energies of 168-basepair-DNA loops for different supercoiling restrains, which are obtained by 3 
iterations as demonstrated in the text. 



ALfc 0.0 0.2 

{E) (in keT) 35.09 36.03 
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In previous literatures, e.g. Refs. p,pl, writhing num- 
ber [Wr) was often used to characterize the tertiary 
structure and handedness of circular DNA, which is de- 
fined as the difference between the linking number Lk 
of two strands and the twisting number Tw of basepairs 
|0. However, Wr has some features which obviously 
hinder itself as the most proper definition of the puck- 
ering degree of DNA loops. For example, for a circular 
DNA wrapping around a sphere, it may be considerably 
puckered, however, Wr definitely equals to zero; when 
DNA takes a small displacement of passage through it- 
self, the puckering degree and aplanarity should keep al- 
most unchanged, however, Wr discontinuously jumps. 

To give a proper definition of the degree of puckering 
of DNA loops, let us at first define a 3 x 3 symmetric 
coordinate tensor: 



Tab — 



3/(?'a(g) - rDa){rb{s) - rob)ds 
/(r(s) -ro)2ds 



a, 6 -1,2, 3, (3) 



where r(s) denotes the axis vector of DNA polymer along 
with its arc-length s, rg — (1/i) / r{s)ds the center of 
mass of the DNA polymer with total arc-length of L, 
ra{b) (s) the projection of the axis vector on a(&)-axis of 3- 
D Cartesian coordinate system. The positive definite ten- 
sor Tab has three eigenvalues T^'s with < Ti < T2 < Ts 
and Ti +T2 + T3 = 3. In fact. Tab represents the config- 
urational ellipsoid of DNA polymer and Tj its i'th major 



axis. So the smallest eigenvalue Ti signifies the degree of 
puckering of the DNA loop from a planar circle and we 
call it "aplanarity of DNA loop" . 

In Fig. 3 is shown the major axes of DNA loop of 
168 bp with different supercoiling stresses, which are 
calculated through the reweighting equation (0) after 
2,000,000 MC runs. To access the handedness of the 
buckling, we also present the values of Wr and ATw 
versus ALfc in Fig. 3d. With supercoiling strain increas- 
ing, the planar circle will become unstable and a con- 
formational transition from circle to figure-8 takes place 
at ALkc ^ —1-2, which manifests itself as an abrupt 
jump in all data of major axes, writhing and twisting 
number, as well as bending and twisting energies (data 
not shown) . This critical value of ALkc is in good agree- 
ment with the analytical results Q, i.e. ALkc = V^^/C, 
where bending persistent length A and twisting persis- 
tent length C are taken as 53 nm and 72.5 nm respec- 
tively in our simulations according to corresponding ex- 
perimental data ||l^. When ALfc continues to increase, 
the figure-8 configuration will become unstable again. 
DNA loop will take diplo- or triple interwound super- 
helical conformation after ALfc is beyond about —2.2 or 
—3.4, which can be most clearly identified in the data of 
aplanarity (see Fig. 3a). However, the latter transitions 
are less abrupt than that of circle to figure-8, as denoted 



by the width of peaks in Fig. 3a. 
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FIG. 3. Averaged major axes as well as writhing and 
twisting number as function of linking difference of 168-base- 
pair-DNA loop, with statistic error bars of MC calculations 
inside the symbols. The lines serve to guide the eye. 

In summary, we showed that short supercoiled DNA 
polymer is a glassy system and canonical Metropolis sim- 
ulation tends to get stuck in some local metastable energy 
basins. A new Monte Carlo algorithm with an explicitly 
analytic weight factor is introduced to solve the problem 
of ergodicity breaking, in which the thermalization is re- 
inforced in both low- and high-energy regions. Compared 
with earlier approach of entropic sampling, the probabil- 
ity weight factor of present algorithm is clearly easier to 
be determined and the implementation is about 2-fold 
CPU time-saving. As a general Monte Carlo method, 
the present approach can be also used in other thermo- 
dynamic systems with frustration. 

To characterize quantitatively the degree of puckering 
of the circular DNA, a new quantity of aplanarity has 
been introduced as the shortest major axis of configura- 
tional tensor of DNA. With the developed Monte Carlo 
method, a quantitative correlation between the buckling 
of DNA loop of 168 basepairs and supercoiling degree is 
presented. Abrupt configurational transition from circle 
to figure-8 takes place at the critical supercoiling stress 



of ALkc (~ 1.2), which is in good agreement with pre- 
vious analytical results. With further increasing linking 
difference, DNA loops will take successively diplo- and 
triple supcrhclical conformation through the decreasingiy 
abrupt transitions. 
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